***************************************************************************************************************************
/*	This do-file
			a. Defines a program that computes a rank preservation test. The steps are as follows:
				1. Compute differences across treatment arms and quartiles. 
				2. The standard errors are computed by randomly assinging synthetic treatment groups, 
				   dividing by quartiles, and computing the synthetic differences. 
				3. A total of J bootstrap replications are done, in order to estimate a distribution of impacts. 
				4. The 5th and the 95th percentiles are taken as lower and upper critical values, such that a 90% 
				   interval is computed. 
			2. The bootstrap intervals are not always equivalent because there is bunching in particular percentiles. 
			   Thus, each time the do-file runs not the exact same observations are included in the same quartiles. The
			   differences, however, are very small. 
*/
***************************************************************************************************************************
// DEFINE PROGRAM: rankbstrp 	
cap program drop rankbstrp
program define rankbstrp
syntax varlist, Reps(integer) [individuals conditional]
// Create Quartiles
	clear matrix	
	set more off
	qui gen group=.
	forvalues i = 0/2 {
		qui sum EL_EGRA_PCA_Index if Study_Arm==`i', d
		qui replace group = 1 if EL_EGRA_PCA_Index<=r(p25) & Study_Arm==`i'
		qui replace group = 2 if (EL_EGRA_PCA_Index>r(p25) & EL_EGRA_PCA_Index<=r(p50)) & Study_Arm==`i'
		qui replace group = 3 if (EL_EGRA_PCA_Index>r(p50) & EL_EGRA_PCA_Index<=r(p75)) & Study_Arm==`i'
		qui replace group = 4 if (EL_EGRA_PCA_Index>r(p75) & EL_EGRA_PCA_Index<.) & Study_Arm==`i'
	}
// COMPUTE MEAN DIFFERENCES FOR ALL COVARIATES (Matrix diff has all the mean differences). 
	mat diff = J(1,8,.)
	local n = 1
	foreach var in `varlist' {
		local c = 1
		forvalues j = 1/4 {
			if "`conditional'" == "" {
				qui reg `var' MT_Program CCT_Program if group==`j'
			}
			else if "`conditional'"!="" {
				qui areg `var' MT_Program CCT_Program if group==`j', a(Group_detailed) 
			}
			mat aux = e(b)
			scalar MT_`var'_`j'= aux[1,1]
			scalar CCT_`var'_`j' = aux[1,2]
		}	
		mat b_`var' = MT_`var'_1,MT_`var'_2,MT_`var'_3,MT_`var'_4,CCT_`var'_1,CCT_`var'_2,CCT_`var'_3,CCT_`var'_4
		mat diff = diff \ b_`var'
	local n = `n'+ 1
	}
	display as result "Differences in Means per Quartile:"
	mat diff = diff[2..`n',1..8]
	mat list diff 
// ESTIMATE BOOSTRAP SE
	// Create Bootstrap sample with replacement
	set seed 651356
	// Start Bootstrap algorithm	
	local matsize = 11000
	set matsize `matsize'	
	local wid = 8 * (`n'-1)
		// Create matrix for exporting
		mat CI_full = J(1,`wid',.)
		foreach var in `varlist' {
			mat CI_full_`var' = J(1,8,.)
		}
	// Algorithm
	tempfile data 
	save `data'
	preserve	
		local rep = 1
		while `rep' <= `reps' {
			use `varlist' EL_EGRA_PCA_Index Group_detailed School_blind Study_Arm using  `data', clear
			// OBTAIN BOOTSTRAP SAMPLE
			display as text "Rank Preservation Bootstrap replication: `rep'"
			qui count
			local obs = r(N)
			qui tab School_blind
			local schools = r(r)
			display _skip(2) "Before sampling:" _newline(1) _skip(5) "Number of individuals = `obs'",  _newline(1) _skip(5) "Number of schools: `schools'"
			// Sample Schools
			bsample, strata(Study_Arm) cluster(School_blind)
			qui count
			local obs = r(N)
			qui tab School_blind
			local schools = r(r)
			display _skip(2) "After sampling schools:" _newline(1) _skip(5) "Number of individuals = `obs'",  _newline(1) _skip(5) "Number of schools: `schools'"
			if "`individuals'" != "" {
				// Sample Individuals
				bsample
				qui count
				local obs = r(N)
				qui tab School_blind
				local schools = r(r)
				display "Sample individuals within schools ON" _newline(1) _skip(2) "After sampling individuals:" _newline(1) _skip(5) "Number of individuals = `obs'",  _newline(1) _skip(5) "Number of schools: `schools'"		
			}
		// Locals for number of observations
			qui count if Study_Arm==0
			local N_CON=r(N)+1
			qui count if Study_Arm==1
			local N_CCT= `N_CON' + r(N)
			qui count if Study_Arm==2
			local N_MT = `N_CCT'+ r(N)-1
			// Sort the data
			qui gen i = runiform()
			sort i
			// Create synthetic Treatments
			qui gen     syn_Study_Arm = 0 in 1/`N_CON'
			qui replace syn_Study_Arm = 1 in `N_CON'/`N_CCT'
			qui replace syn_Study_Arm = 2 in `N_CCT'/`N_MT'
			// Create synthetic treatment status
			qui gen MT_Program = (syn_Study_Arm==1)
			qui gen CCT_Program = (syn_Study_Arm==2)
			// Define Quartile values 
			qui gen group=.
			forvalues i = 0/2 {
				qui sum EL_EGRA_PCA_Index if syn_Study_Arm==`i', d
				qui replace group = 1 if EL_EGRA_PCA_Index<=r(p25) & syn_Study_Arm==`i'
				qui replace group = 2 if (EL_EGRA_PCA_Index>r(p25) & EL_EGRA_PCA_Index<=r(p50)) & syn_Study_Arm==`i'
				qui replace group = 3 if (EL_EGRA_PCA_Index>r(p50) & EL_EGRA_PCA_Index<=r(p75)) & syn_Study_Arm==`i'
				qui replace group = 4 if (EL_EGRA_PCA_Index>r(p75) & EL_EGRA_PCA_Index<.) & syn_Study_Arm==`i'
			}
			*tab group // To check differences in bunching.
			// Take differences
			mat CI_full_`rep' = J(1,1,.)
			foreach var in `varlist' {
				forvalues j = 1/4 {
					if "`conditional'" == "" {
						qui reg `var' MT_Program CCT_Program if group==`j'
					}
					else if "`conditional'"!="" {
						qui areg `var' MT_Program CCT_Program if group==`j', a(Group_detailed) 
					}
					mat aux = e(b)
					mat `var'_MT_CI_`j'= aux[1,1]
					mat `var'_CCT_CI_`j' = aux[1,2]
				}
				// For full replications	
				mat `var' = `var'_MT_CI_1,`var'_MT_CI_2,`var'_MT_CI_3,`var'_MT_CI_4,`var'_CCT_CI_1,`var'_CCT_CI_2,`var'_CCT_CI_3,`var'_CCT_CI_4	
				mat coln `var' = "MT_`var'_1" "MT_`var'_2" "MT_`var'_3" "MT_`var'_4" "CCT_`var'_1" "CCT_`var'_2" "CCT_`var'_3" "CCT_`var'_4" 
				mat CI_full_`rep' = CI_full_`rep',`var'	
				// Extract matrix with values for each variable
				mat CI_full_`var' = CI_full_`var' \ `var'
			}
			// Extract Matrix with all bootstrap replications for all covariates	
			local wid2 = `wid' + 1	
			mat CI_full_`rep' = CI_full_`rep'[1,2..`wid2']	
			mat CI_full = CI_full \ CI_full_`rep'
			mat drop CI_full_`rep'
			local rep = `rep' + 1
		}
		// Extract Values
		mat CI = J(1,16,.)
		mat pval = J(1,8,.)
		foreach var in `varlist' {
			clear
			qui svmat CI_full_`var'
			qui drop in 1
			forvalues t = 1 / 8 {
				_pctile CI_full_`var'`t',nq(100)
				scalar a = round(r(r5),0.001)
				scalar b = round(r(r95),0.001)
				mat aux_`t' = a, b
				mat coln aux_`t' = "Low" "High"
				gen above_est = abs(CI_full_`var'`t')>abs(b_`var'[1,`t'])
				sum above_est
				scalar pval_`t' = r(mean)
				drop above_est				
			}
			mat CI_`var' = aux_1,aux_2,aux_3,aux_4,aux_5,aux_6,aux_7,aux_8
			mat CI = CI \ CI_`var'
			mat pval_`var' = pval_1,pval_2,pval_3,pval_4,pval_5,pval_6,pval_7,pval_8
			mat pval = pval \ pval_`var'
		}
		mat CI = CI[2..`n',1..16]
		mat pval = pval[2..`n',1..8]
		mat total = diff, CI, pval
	restore
// Create New Data Set
	local n = 1
	qui gen names = ""
	foreach var in `varlist' {
		local lab: variable label `var'
		qui replace names = "`lab'" in `n'
		local n = `n'+ 1
	}
	drop if names==""
	mat coln total = "MT_p25"           "MT_p25_p50"       "MT_p50_p75"       "MT_p75" ///
					 "CCT_p25"          "CCT_p25_p50"      "CCT_p50_p75"      "CCT_p75" ///
					 "lci_MT_p25"       "hci_MT_p25"       "lci_MT_p25_p50"   "hci_MT_p25_p50" ///
					 "lci_MT_p50_p75"   "hci_MT_p50_p75"   "lci_MT_p75"       "hci_MT_p75" ///
					 "lci_CCT_p25"      "hci_CCT_p25"      "lci_CCT_p25_p50"  "hci_CCT_p25_p50" ///
					 "lci_CCT_p50_p75"  "hci_CCT_p50_p75"  "lci_CCT_p75"      "hci_CCT_p75" ///
					 "pval_MT_p25"      "pval_MT_p25_p50"  "pval_MT_p50_p75"  "pval_MT_p75" ///
					 "pval_CCT_p25"     "pval_CCT_p25_p50" "pval_CCT_p50_p75" "pval_CCT_p75" 
	qui keep names
	qui svmat total, names(col)
	// Put Stars
	local perc MT_p25 MT_p25_p50 MT_p50_p75 MT_p75 CCT_p25 CCT_p25_p50 CCT_p50_p75 CCT_p75
	foreach var of local perc {
		qui gen star_`var' = "*" if pval_`var'< 0.1
		qui replace star_`var' = "**" if pval_`var'< 0.05
		qui replace star_`var' = "***" if pval_`var'< 0.01
		*qui gen star_`var' = "*" if `var'< lci_`var' | `var'>hci_`var'
	}
	// I will append the CI under the point estimates. Therefore, the names 
	// of the variables have to be changed. 
	preserve
		keep names lci_* hci_*
		gen order = _n
		foreach var of local perc {
			rename lci_`var' `var'
			rename hci_`var' h_`var'
		}
		tempfile aux
		save `aux'
	restore
	// Change names of point estimates
	qui drop lci_* hci_*
	qui gen order = _n
	foreach var of local perc {
		qui gen h_`var' = . 
	}
	// Append point estimates over the standard errors. 
	qui append using `aux', gen(ci)
	// Add stars and fix data
	qui gen aux_1 = "[" if ci==1
	qui gen aux_2 = "]" if ci==1
	foreach var of local perc { 
		qui tostring `var', replace force format(%9.3f)
		qui tostring h_`var', replace force format(%9.3f)
		qui replace `var' = `var'+star_`var' if ci==0
		qui replace `var' = aux_1 + `var'+ ";"+h_`var' + aux_2 if ci==1
	}
	sort order ci
	qui replace names = "" if ci==1
	qui drop ci	star_* h_* aux_* order pval*
	// Labels
	label var MT_p25       "0-25th Percentile"
	label var MT_p25_p50   "25-50th Percentile"
	label var MT_p50_p75   "50-75th Percentile"
	label var MT_p75       "75-100th Percentile"
	label var CCT_p25      "0-25th Percentile"
	label var CCT_p25_p50  "25-50th Percentile"
	label var CCT_p50_p75  "50-75th Percentile"
	label var CCT_p75      "75-100th Percentile"	
end